Ising model on two connected Barabasi- Albert networks 



O 

o 



Krzysztof SucheckQ and Janusz A. HolyslQ 
Faculty of Physics and Center of Excellence for Complex Systems Research 
Warsaw University of Technology 
Koszykowa 75, PL-00-662 Warsaw, Poland 
(Dated: February 6, 2008) 

We investigate analytically the behavior of Ising model on two connected Barabasi- Albert net- 
works. Depending on relative ordering of both networks there are two possible phases corresponding 
to parallel or antiparallel alingment of spins in both networks. A difference between critical tem- 
peratures of both phases disappears in the limit of vanishing inter-network coupling for identical 
networks. The analytic predictions are confirmed by numerical simulations. 



C/3 



I 

o 
o 



> 
m 
o\ 
\o 
cn 
o 

o 



I 

o 
o 



X 



PACS numbers: 05.50.-|-q, 



i.75.-k, 



i.75.Fb 



I. INTRODUCTION 

Isinffmodel on Barabasi- Albert (B-A) scale-free net- 
work |l| has been investigated both numerically as 
well as analytically Q and it has been shown that the 
critical temperature of such a system is proportional to 
the logarithm of a total spins number N. Similar Ising 
models have been investigated for a general class of ran- 
dom scale-free graphs ^ ,5, 6] and it has been shown 
that critical temperatures of such systems depend sub- 
stantially on a characteristic exponent 7 describing the 
probability distribution of node degrees. Other studies 
of this model include investigations of antiferromagnetic 
interactions 0, dynamics on directed networks ^ and 
critical properties of spin-glass . 

Besides a large interest in physical properties of Ising- 
like models they seem also to be important for opinion 
formation modeling 0, 0, 0, ^|. The Ising model 
exhibits a majority rule dynamics - the feature that of- 
ten can be found in social systems, where a given per- 
son changes his/her opinion to fit to a majority of his 
neighbors. Since it is common that social networks have 
modular structure of weakly coupled clusters it is of 
particular interest to see how the Ising model behaves in 
the case of two interacting complex networks. While ge- 
ometrical properties of interconnected complex networks 
have been studied before [l5l |. the dynamics in such sys- 
tems has not been explored throughly. 
We start with analytical investigations of Ising model on 
two interconnected Barabasi- Albert networks, and then 
show results of numerical simulations confirming our an- 
alytical studies. 



II. MODEL 

Our model considers Ising spins on a B-A network. 
Each node of the network has a spin. We study only 



ferromagnetic interactions existing between directly con- 
nected spins. 




B 



FIG. 1: Two connected B-A networks. A few nodes from each 
network are shown. The intra- network degrees kAA and ksB 
as well as inter-network degrees kAB and fcss for two sample 
nodes are presented. 

The B-A model is a model of a growing network Q]. 
One starts with m fully connected nodes, and adds new 
nodes to the network. Each new node creates m connec- 
tions to the extisting network. The probability that a 
connection will be made to a node i is proportional to 
its degree ki. This results in a scale- free network, with a 
degree distribution P{k) ^ k~~^. 

We assume that two B-A networks are connected by Eab 
links (FigQ . Each of these links connects a node in net- 
work A with a node in network B. The chance to choose 
a given node as the end of the inter-network link is pro- 
portional to the intra-nctwork node degree. This means 
that for a small number of links Eab, an inter-network 
degree kAB of a node in the network A is proportional 
to a intra-network node degree kAA- A similar relation 
holds for degrees in the network B. 



III. ANALYTIC CALCULATIONS 
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The problem of the Ising model in a single B-A net- 
work was already solved analytically by Bianconi by an 
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appropriately tailored mean-field approach ■ We use a 
similar approach for the problem of two connected net- 
works. 

The Hamiltonian of the Ising model for a single B-A net- 
work can be written as 
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where Si, Sj — ±1 are spins of nodes a constant Jij 
is a ferromagnetic coupling between them and hi is an 
external field acting on spin i. The coupling constants 
Jij equal to a positive constant if spins are connected, 
and are zero otherwise. 

The exact solution for the average spin {si) in a single 
network can be written as 



(si) = ( tanh I (3 ^ JijSj 




(2) 



where /3 = l/T, the temperature T is measured in units 
of inverse Boltzmann constant l/fcs and averaging is over 
the canonical ensemble. If we now consider an aver- 
age over all possible realizations of B-A networks, then 
(Jij) = Jkikj/E. We use the mean field approximation, 
taking {Jij) in place of Jij in our equation. Since i and j 
are ordered pairs, the total number of pairs E is twice the 
number of edges in the network. If we take the external 
field hi equal to zero, the Eq. ^has the following form 



{si) = i tanh 



(3) 



Now we consider a pair of coupled networks A and 
B. The parameters describing both networks can be 
split into four groups - two describe internal properties 
of each network, and two describe network-network 
interactions. We introduce the following notation: sai 
and SBi are spins in networks A and B, Jaa, Jbb are 
coupling constants between spins in networks A and 
B respectively. Jab = Jba arc the coupling constants 
between spins in different networks, kAAi and ksBi 
are intra-network node degrees, kABi and ksAi are 
inter- network node degrees, Eaa and Ebb are twice the 
total numbers of all intra-network links in A and B, 
Eab = Eba is the number of links between the networks. 

Now we extend the Eq. O introducing the influence of 
the second network. This way we obtain two equations 
for average spins in every network. Since we are inter- 
ested in critical properties, where average spins are close 
to zero, we can approximate the hyperbolic tangent by a 
linear function and use a standard mean-field approach. 
As result we get 



{sAi) = l3JAAkAAi 
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kAAj jsAj) 

Eaa 



To get a relation for the system critical temperature we 
need to have a self-consistent equations for order param- 
eter. The case of a single network required introduction 
of only single weighted spin 5" = J2i ki {si) /E, where the 
(si) is mean- field average for a given spin i. In the case 
of two connected networks, we need to consider four such 
weighted spins Saa, Sbb, Sab and Sba- 



Saa = ^ kAAi {sAi) I Eaa 

i 

Sbb — ^ kBBi {sBi) / Ebb 

i 

Sab — kABi (sAi) / Eab 

i 

Sba = kBAi {sBi) /Eba 



(6) 
(7) 
(8) 
(9) 



Saa and Sb b hold the same meaning as for a single net- 
work. Sab is the mean weighted spin of the network A 
observed by spins in the network B, while Sba is a mean 
weighted spin of the network B observed by spins in the 
network A. 

Eas. HlUll3|) were received from Eqs . (|4l5f) by multiplying 
by appropriate factors (see Eas lt)l9|l and summing over i. 
These four equations contain only four weighted spins as 
unknown collective variables and are approximate mean- 
field description of the system close to a critical point. It 
follows we receive 



Saa ~ PJaaSaa -I- PJbaSba 

I i 

Sbb ~ PJbbSbb ^ ' J?^' + PJabSab ^ ^ 



Sab ~ PJaaSaa 
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Sba = PJbbSbb ^ ^ 
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(10) 



(11) 



+ PJBASBAy ^{12) 

^-^ Eab 



+ /3JAflSAfl^^(13) 



If we assume that kABi = PAkAAi and kBAi = PBkBBi, 
what means that the number of links outside the network 
is proportional to the number of links within the network, 
we can greatly simplify our four equations. This can be 
done when one takes into account the way we create inter- 
network links in our model. 

Using this assumption, we do not need to consider the 
cross-network weighted spins Sab and Sba as they are 
proportional to Saa and Sbb- Now our first two equa- 
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tions become 
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(15) 



where 5*^ = Saa and Sb = Sbb- The equation array 
can be written as a single matrix equation. 

S = (3AS (16) 



where S is a vector 



Sa 
Sb 



describing the state of the 



system and A is a matrix describing effective interaction 
strengths between spins belonging to the same or to dif- 
ferent networks 



Aaa Asa 
Aas ^bb 



PaJab\^ 



PbJbaYe_ 



7J 



(17) 



In the case of a single network A, solutions other than 



Sa — can exist only if PJaa 



(kAA) 



> 1 In the case 



of two coupled networks, this condition corresponds to 
an eigenvalue of Eq^l greater than 1 . The eigenvalues 
are 



/3 



(^Aaa + ^bb ± V (Aaa - Ags f + 4AflAAAs) (18) 

Comparing these eigenvalues with 1, we get the follow- 
ing critical temperatures 



^ Aaa + Abb ± \/ {Aaa - Ass)^ + 4Ai3AAAi3 

Tc± = 2 (.^^^ 

Since the diagonal elements of A, Aaa and Abb are 
critical temperatures Tca,Tcb for separate networks we 
can write the critical temperatures for the coupled system 
as 



T,A + T,b ± V {TcA - T,b) + ^AbaAab 
Tc±^ ^- ^ (20) 

To better understand the meaning of these solutions, 
we introduce the following variables 



A^{T,a + T,b)/2 
V = {T,A - T,b)I'2 
C = \/ AbaAab 



(21) 
(22) 
(23) 



The value A ("average") describes an average critical 
temperatures of the networks, T) ("difference") is the dif- 
ference between critical temperatures of both networks, 
C ("coupling") describes a strength of internetwork 
interactions. 

Using this notation the critical temperatures can be 
written shortly as 



(24) 



Let us now consider eigenvectors associated with X±. 
They are proportional to the magnetization of both net- 
works that appears below a given critical temperature 
and disappears above it. The unnormalized eigenvectors 
are 



S+ = 



1 



Aba 



(25) 



The eigenvector S_ has opposite signs of its components 
and corresponds to networks ordered with antiparallel 
weighted spins, while the eigenvector S+ has the same 
signs of the components and corresponds to networks or- 
dered with parallel weighted spins. 

In the limit of vanishing internetwork coupling (C — 0) 
the eigenvalues are simply the diagonal elements of the 
matrix Xa — Tca, Xb — Tcb and the associated nor- 



malized eigenvectors are — 



Sb - 



This 



means that in this limit two stable states of the system 
correspond to the ordering of just one of the networks, 
and there is no relation between the order in each net- 
works. It shows our approach gives correct results in this 
specific case. 

Now let us investigate the inequality conditions for ex- 
istence of solutions S± . From the condition A± > 1 we 
have 



A-T> yJV^ + C2 



(26) 
(27) 




FIG. 2: The area of parameters C and "D where the ordered 
states S± can exist. The order appears in the dashed areas. 
The state S- of antiparallel ordering corresponds to the left 
picture, the state S+ or parallel ordering corresponds to the 
right picture. The radius of the circle is the difference T — A. 
Note that solution S_ does not exist for A <T and solution 
S+ exist for any C and T) for A > T. In the case C = the 
networks do not interact at all, and the solutions lose their 
normal meaning. 

The meaning of these inequalities is presented at 
FigH 

If we consider networks of the same size, the dependence 
of critical temperatures on inter-network interaction 
strength C is linear, and both critical temperatures are 
the same for C ^ 0. In such case the system critical 
temperature = A = T^a = Tcb (Fig'EIl- 

The analytic results hold true for any random network, 
where the probability of a link existing between any two 
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nodes i and j is proportional to the product kikj. This 
is the only assumption about network structure we have 
used, so any networks where the condition is fulfilled (i.e. 
random networks) is described by our analysis. 
Our numeric calculations found in the following section, 
correspond to the specific case of B-A network and con- 
stant coupling Jaa = Jab = Jba = Jbb = J ■ We can 
write the critical temperatures as follows 

r.± = 2 ^^"^ )' + P-^P^T.^'^^^ (28) 

where T^a = J(™a/2) IniV^ and Tcb = Jims/'i) InNs- 
The values of and ps are not independent, and are 
connected with the number of links between networks 

PaEaa = PbEbb = Eab- 



correlated. We find the weighted spin S = Sa + Sb for 
each temperature T and average it over 100 network 
realizations. Simulations are performed for different 
numbers of inter-network finks Eab = 1000, 2000, 3000 
and 4000 that were attached preferentially to fulfill the 
assumptions of our model (see the discussion at the end 
of SectHH). 

For low temperatures T the system is ordered. As 
the temperature increases the average weighted spin 
S decreases. When the temperatures increase over 
Tc+, the ferromagnetic ordered state of both networks 
disappears, i.e. both networks become paramagnetic 
(Fig©. 
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FIG. 3: The dependence of temperatures Tc± on number of 
links between networks Eab- The dashed lines are analytic 
predictions fEa ll9^ . The symbols are numeric results. Circles 
correspond to Tc- . Squares and diamonds correspond to Tc+ 
and are calculated from susceptibility x ~ {Sh — So)/h. Trian- 
gles correspond to Tc+ and are calculated from susceptibility 
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FIG. 4: The dependence of the average weighted spin S = 
Sa + Sb and its susceptibility for small external field h = 
0.05J in the case of two B-A networks with Eab = 1000 
connections between them. An initial condition for each 
temperature is a completely ordered system. The full sym- 
bols depict S, the X symbols correspond to susceptibility 
X — {So.05 ~ So) /h (the lines are just to guide eye), the 
empty symbols are 30-point running average of the suscep- 
tibility. The parabolic fit was used to find the susceptibility 
maximum. 



IV. NUMERIC RESULTS 

Our analytic calculations show the existence of two 
different ordered states and estimate values of two 
critical temperatures where these states disappear. 
Below we investigate numerically a case of two cou- 
pled B-A networks with the same number of nodes 
Na = Nb = 5000 and finks Eaa = Ebb = 40000 
{ruA ^TUB = 5, (kAA) = {kBs) = 10). 

We run Ising dynamics on these networks, setting the 
following initial condition: all spins in both networks 
have the same value SAi = ssi = +1. We allow the 
system to relax for r = 20 time steps, perform averaging 
for T = 20 time steps, then increase the temperature 
and start from the same initial condition as before. 
This way, results for different temperatures are not 



Finding the exact value of critical temperature Tc+ 
from numerical simulations is not straightforward. If one 
observes the dependence of the weighted magnetization 
on rising temperature and tries to fit the magnetization 
decay to a linear or to an exponential function the results 
strongly depend on relaxation time r. To overcome this 
problem we observed the temperature dependence of 
the system susceptibility x- In fact, by comparison to 
standard models of magnetic systems, one can expect 
that the initial susceptibility diverges at T = Tc+. In our 
finite system we are looking simply for the maximum of 
X- To estimate x we are using two methods. First we 
compare average weighted spin S for a small external 
field h = 0.05J and the value of 5*0 with no external 
field. It follows X = {Sh — So) /h. Because such results 
are strongly fiuctuating as a function of system history 
and temperature (Fig^J, we calculate running average 
over 30 temperature points and find the maximum of x 
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by fitting a parabolic curve. The top of the parabola 
corresponds to the position of the critical temperature 
Tc+- We found that these values are independent on the 
relaxation time r used in our numerical experiment. The 
second method of finding the critical temperature Tc+ is 
observation of the time average {S"^) - where we 

average over one relaxation period r. The magnitude 
of the fluctuations is proportional to the suceptibility 
X ~ ('5'^) — {S)^ according to the fluctuation-dissipation 
theorem IQj. Similarly to previous method, we cal- 
culate running average over 10 points and find the 
maximum. The values are shifted by a constant value 
comparing to analytic results and do not fluctuate as 
much as those obtained from the first method (see Fig|3J|. 

Now we consider the same networks with the following 
initial condition for each temperature: spins in both net- 
works are ordered antiparallel SAi = = —1. The 
relaxation algorithm and the measurement procedure is 
the same as for parallel spin case, however now we con- 
sider absolute values of weighted spin of the whole sys- 
tem \S\ = \Sa + Sb\- For low temperatures T both net- 
works remain ordered in opposite directions (Fig[SJ and 
as result an average weighted spin value |S'| = jS'yi -I- ^sj 
fluctuates around 0. When the temperature increases 
over Tc- , the state of antiparallely ordered networks dis- 
appears and networks start to order in a parallel fash- 
ion. An example of such a scenario is presented at FigEl 
When temperature further increases and T > Tc-^- the 
system becomes paramagnetic. 



declining part of the curve S{T) (linear fit was used 
for the rising part and exponential fit for the declining 
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FIG. 5: The temperature dependence of the average weighted 
spin value |S| = \Sa + Sb\ in the case of two B-A networks 
with Eab = 4000 connections between them. Initial con- 
ditions for each temperature are two fully ordered networks 
with opposite spins. The symbols correspond to numerical 
data and grey lines are approximations of magnetization be- 
havior below and above critical temperature. 

We found numerically the critical temperature Tc_ 
from the intersection of extrapolations of rising and 
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FIG. 6: An example of time evolution of magnetization for 
two networks of size A''^ = Nb ~ 5000 and (k) = 4 with 
Eab = 500 links between them in temperature T = 30 (scaled 
in kB/J)- At this temperature, the antiparallel ordered state 
is not stable, and the initial condition of networks ordered in 
such fashion evolves towards a stable parallel ordered state. 



part) (FigO. This point slightly depends on relaxation 
times T. For large r the effect of finite system size 
can be easily observed and the system jumps from an 
antiparallel state to a parallel one that has a lower energy. 

The results for dependence of Tc_ and Tc+ on the num- 
ber of links between networks Eab agree with the ana- 
lytic calculations (Fig|3Jl. 



V. CONCLUSIONS 

In the system of two B-A network, the Ising model 
possesses two low-temperature stable states — both net- 
works ordered parallel or antiparallel. It follows there are 
two critical temperatures corresponding to the disappear- 
ance of these two stable states — Tc_ and Tc+ . They are 
placed symmetrically around the average of critical tem- 
peratures of separate networks. The difference between 
them depends on density of inter-network links and the 
difference between critical temperatures of separate net- 
works. The analytic calculations agree with performed 
numeric simulations. 
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